# Paper: Is it still the economy? Economic voting in polarized politics
# County-level model
# Replication: Figure 4
remove(list = ls())

library(conflicted)
library(tidyverse)
conflict_prefer("filter","dplyr")
library(ggplot2)
library(dplyr)
library(here)
conflict_prefer("here", "here")
library(systemfit)

load(here("Final_Paper", "R&R", "CountyLevel.RData"))


lia1 <- inc_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff +
  inc_abs_lag + state + rep_inc + incumbent_cand
lia2 <- opp_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) +
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff +
  opp_abs_lag + state + rep_inc + incumbent_cand
Systemia <- list(lia1, lia2)
modelia <- systemfit(Systemia, method = "SUR", data = model_data)
summary(modelia)

lio1 <- inc_opp ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff + 
  inc_opp_lag + state + rep_inc + incumbent_cand
lio2 <- abs_opp ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff + 
  abs_opp_lag + state + rep_inc + incumbent_cand
Systemio <- list(lio1, lio2)
modelio <- systemfit(Systemio, method = "SUR", data = model_data)
summary(modelio)


# abstention baseline
(estia <- cbind(Estimate = coef(modelia), confint(modelia, level = 0.95)))
estia
length(estia[,1])
summary(modelia)
estia[2,]
eq1iau <- estia[2, ] #unemployment rate
eq1iau
estia[74, ] #148/2
eq1iaup <- estia[74, ] #interaction
eq1iaup
eq1ia <- rbind(eq1iau, eq1iaup)
eq1ia

estia[76, ]
summary(modelia)
eq2oau <- estia[76, ] #unemployment rate
eq2oau
estia[148, ]
eq2oaup <- estia[148, ] #interaction
eq2oaup
eq2oa <- rbind(eq2oau, eq2oaup)
eq2oa

variables <- c("Unemp.Rate", "Unemp.Rate*Polarization")

eq1ia_data <- data.frame(variables, eq1ia) %>% 
  mutate(model = "Equation 1") %>% 
  mutate(equation = "Inc/Abs")
glimpse(eq1ia_data)
eq2oa_data <- data.frame(variables, eq2oa) %>% 
  mutate(model = "Equation 2") %>% 
  mutate(equation = "Opp/Abs")
glimpse(eq2oa_data)


# opposition baseline 
(estio <- cbind(Estimate = coef(modelio), confint(modelio, level = 0.95)))
estio
summary(modelio)
estio[2, ]
eq1iou <- estio[2, ] #unemployment rate
eq1iou
estio[74, ]
eq1ioup <- estio[74, ] #interaction
eq1io <- rbind(eq1iou, eq1ioup)
eq1io

eq1io_data <- data.frame(variables, eq1io) %>% 
  mutate(model = "Equation 3") %>% 
  mutate(equation = "Inc/Opp")

model_datai <- rbind(eq1ia_data, eq2oa_data, eq1io_data)
glimpse(model_datai)

# standard errors
summary(modelia)
coef(summary(modelia))[2,2]
coef(summary(modelia))[74,2]    
se_iau <- coef(summary(modelia))[2,2]
se_iaup <- coef(summary(modelia))[74,2]
se_oau <- coef(summary(modelia))[76,2]
se_oaup <- coef(summary(modelia))[148,2]
se_iou <- coef(summary(modelio))[2,2]
se_ioup <- coef(summary(modelio))[74,2]

vcov_data <- as.data.frame(vcov(modelia))
vcov(modelia)
sqrt(vcov(modelia)[2, 2])
coef(summary(modelia))[2,2]
vcov(modelia)[2, 74]
coviaup <- vcov(modelia)[2, 74]
sqrt(vcov(modelia)[76, 76])
coef(summary(modelia))[76,2]
covoaup <- vcov(modelia)[76, 148]
covioup <- vcov(modelio)[2, 74]

remove(vcov_data)
se_i <- rbind(se_iau, se_iaup, se_oau, se_oaup, se_iou, se_ioup)
model_datai <- cbind(model_datai, se_i)
glimpse(model_datai)
sqrt(vcov(modelia)[2, 2])
sqrt(vcov(modelia)[74, 74])

model_datai <- model_datai %>% 
  mutate(covb = ifelse(variables == "Unemp.Rate*Polarization" & 
                         equation == "Inc/Abs", coviaup, 0),
         covb = ifelse(variables == "Unemp.Rate*Polarization" & 
                         equation == "Opp/Abs", covoaup, covb),
         covb = ifelse(variables == "Unemp.Rate*Polarization" & 
                         equation == "Inc/Opp", covioup, covb))
glimpse(model_datai)

# levels of polarization (0.6, 0.7, 0.85, 0.9)
model_datai06 <- model_datai %>% 
  mutate(Estimate = ifelse(variables == "Unemp.Rate*Polarization",
                           0.6*Estimate, Estimate)) %>% 
  mutate(polarization = 0.6) 
model_datai08 <- model_datai %>% 
  mutate(Estimate = ifelse(variables == "Unemp.Rate*Polarization",
                           0.85*Estimate, Estimate)) %>% 
  mutate(polarization = 0.85)
model_datai <- model_datai %>% 
  mutate(polarization = 1)

model_datai <- rbind(model_datai, model_datai06, model_datai08)
glimpse(model_datai)
table(model_datai$polarization)

# variables for variance of marginal effect
model_datai <- model_datai %>% 
  mutate(varb = ifelse(variables == "Unemp.Rate",
                       se_i^2, polarization^2*se_i^2)) 
glimpse(model_datai)

model_datai2 <- model_datai %>% 
  mutate(variables = ifelse(variables == "Unemp.Rate", "Unemp.Rate*Polarization", variables)) %>% 
  group_by(equation, polarization, variables) %>% 
  summarise(Estimate = sum(Estimate),
            vars = sum(varb),
            covb = sum(covb)) %>% 
  ungroup()
glimpse(model_datai2)

# standard error of marginal effects
model_datai2 <- model_datai2 %>% 
  mutate(se = sqrt(vars + 2*polarization*covb))
model_datai2 <- model_datai2 %>% 
  select(equation:Estimate, se)

glimpse(model_datai)
model_datai1 <- model_datai %>% 
  filter(polarization == 1) %>% 
  filter(variables == "Unemp.Rate") %>% 
  mutate(polarization = 0)
model_datai1 <- model_datai1 %>% 
  select(equation, variables, polarization, Estimate, se = se_i)

model_datai <- rbind(model_datai1, model_datai2)
glimpse(model_datai)
model_datai <- model_datai %>% 
  mutate(icii = Estimate+(1.96*se)) %>% 
  mutate(ici = Estimate-(1.96*se))

glimpse(model_datai)

model_datai <- model_datai %>%
  mutate(polarization = as.factor(polarization)) %>%  
  mutate(polarization = factor(polarization, 
                               levels = c("0",
                                          "0.6",
                                          "0.85",
                                          "1")))
model_datai <- model_datai %>%
  mutate(equation = factor(equation, 
                           levels = c("Inc/Opp", "Inc/Abs", "Opp/Abs"))) %>% 
  mutate(equation = ifelse(equation == "Inc/Opp", "log(Incumbent/Opposition)", equation),
         equation = ifelse(equation == "2", "log(Incumbent/Abstention)", equation),
         equation = ifelse(equation == "3", "log(Opposition/Abstention)", equation)) %>% 
  mutate(equation = factor(equation, 
                           levels = c("log(Incumbent/Opposition)", 
                                      "log(Incumbent/Abstention)", 
                                      "log(Opposition/Abstention)")))
glimpse(model_datai)
table(model_datai$polarization)

fig_incopp <- model_datai %>% 
  filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  filter(equation == "log(Incumbent/Opposition)") %>% 
  filter(polarization == "0.6" | polarization == "0.85") %>% 
  ggplot(aes(color = polarization, shape = polarization, fill = polarization)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = .75) +
  geom_pointrange(aes(x = equation, y = Estimate, ymin = ici,
                      ymax = icii),
                  lwd = .8, position = position_dodge(width = .5)) + 
  theme_bw() + ylim(-0.05, 0) +
  ylab("Marginal Effect") + xlab(NULL) + #ylim(-0.1, 0.05) +
  scale_color_manual(name = "Polarization",
                     values = c("#009999", "#990000"), 
                     labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization",
                     values = c(24, 25), 
                     labels = c("0.60", "0.85")) +
  scale_fill_manual(name = "Polarization",
                    values = c("#009999", "#990000"), 
                    labels = c("0.60", "0.85")) +
  ggtitle("(a) log(Incumbent/Opposition)") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 13),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 13))
fig_incopp

fig_incabs <- model_datai %>% 
  filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  filter(equation == "log(Incumbent/Abstention)") %>% 
  filter(polarization == "0.6" | polarization == "0.85") %>% 
  ggplot(aes(color = polarization, shape = polarization, fill = polarization)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = .75) +
  geom_pointrange(aes(x = equation, y = Estimate, ymin = ici,
                      ymax = icii),
                  lwd = .8, position = position_dodge(width = .5)) + 
  theme_bw() + #ylim(-0.05, 0) + 
  scale_y_continuous(limits = c(-0.05, 0)) +
  ylab(NULL) + xlab(NULL) + #ylim(-0.1, 0.05) +
  scale_color_manual(name = "Polarization",
                     values = c("#009999", "#990000"), 
                     labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization",
                     values = c(24, 25), 
                     labels = c("0.60", "0.85")) +
  scale_fill_manual(name = "Polarization",
                    values = c("#009999", "#990000"), 
                    labels = c("0.60", "0.85")) +
  ggtitle("(b) log(Incumbent/Abstention)") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 13),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 13))
fig_incabs

fig_oppabs <- model_datai %>% 
  filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  filter(equation == "log(Opposition/Abstention)") %>% 
  filter(polarization == "0.6" | polarization == "0.85") %>% 
  ggplot(aes(color = polarization, shape = polarization, fill = polarization)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = .75) +
  geom_pointrange(aes(x = equation, y = Estimate, ymin = ici,
                      ymax = icii),
                  lwd = .8, position = position_dodge(width = .5)) + 
  theme_bw() + 
  scale_y_continuous(limits = c(-0.012, 0.05)) +
  ylab("Marginal Effects") + xlab(NULL) + #ylim(-0.1, 0.05) +
  scale_color_manual(name = "Polarization",
                     values = c("#009999", "#990000"), 
                     labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization",
                     values = c(24, 25), 
                     labels = c("0.60", "0.85")) +
  scale_fill_manual(name = "Polarization",
                    values = c("#009999", "#990000"), 
                    labels = c("0.60", "0.85")) +
  ggtitle("(c) log(Opposition/Abstention)") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 13),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 13))
fig_oppabs

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
multiplot(fig_incopp, fig_incabs, cols = 2)







